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SECTION  I 
INTRODUCTION 


1.1  OVERVIEW  AND  TECHNICAL  APPROACH 

The  extended  Kalman  filter  (EKF)  is  a general  non- 
linear estimation  technique  which  can  be  applied  to  the  prob- 
lem of  estimating  the  aerodynamic  coefficients  of  a body  from 
free-flight  measurements  of  its  motion.  The  basic  technique 
can  handle  a wide  range  of  nonlinear  aerodynamic  models  and 
trajectory  measurement  systems,  and  it  accounts  for  random 
errors  in  the  trajectory  measurements,  stochastic  disturbances 
to  the  body's  motion,  and  a priori  information  about  the  aero- 
dynamic coefficients  to  be  estimated.  This  report  describes 
an  application  of  the  EKF  to  the  estimation  of  the  aerodynamic 
coefficients  of  a six-degree-of-freedom  rotationally  symmetric 
rigid  body,  based  on  discrete  free-flight  trajectory  measure- 
ments, such  as  those  made  at  the  Air  Force  Armament  Laboratory 
(AFATL)  Aerobal 1 ist i c Research  Facility  test  range.  These 
measurements  consist  of  three  spatial  positions  and  three 
angular  orientations,  relative  to  a fixed  inertial  coordinate 
system,  and  time-of-f 1 ight  at  50  downrange  positions  along 
the  trajectory.  The  algorithm  incorporates  a stochastic  mea- 
surement model  that  approximates  the  conditions  which  exist 
at  the  test  range. 

In  addition  to  estimates  of  the  projectile  aero- 
dynamic coefficients,  the  EKF  algorithm  provides  an  estimate 
of  the  rms  error  associated  with  each  parameter  estimate. 
Algorithm  performance  is  evaluated  by  estimating  the  aero- 
dynamic coefficients  from  synthetic,  computer-generated  tra- 
jectory data.  These  data  are  derived  from  the  same  projectile 
dynamic  model  that  is  used  to  design  the  filter  algorithm. 

This  report  presents  evaluation  results  including  assessments 
of  the  accuracy  of  the  estimates;  the  consistency  between  the 
estimation  errors  and  their  standard  deviations  computed  by 
the  filter;  and  the  sensitivity  to  projectile  trajectory, 
measurement  noise  level,  and  initial  conditions. 


1.2  ORGANIZATION  OF  THE  REPORT 

Section  II  of  this  report  summarizes  and  discusses 
the  equations  of  the  extended  Kalman  filter  and  provides  the 
mathematical  basis  for  its  application  to  the  estimation  of 
aerodynamic  coefficients.  Section  III  describes  the  details 
of  the  projectile  dynamic  model  and  illustrates  the  applica- 
tion of  the  EKF  technique.  Performance  and  sensitivity  study 
results  are  given  in  Section  IV,  and  Section  V summarizes  the 
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conclusions  of  the  study  and  provides  suggestions  for  future 
related  investigations  and  applications  of  the  EKF  parameter 
estimation  technique.  Section  VI  presents  the  recommendations 
resulting  from  this  effort.  Finally,  Appendix  A gives  some 
of  the  detailed  equations  required  to  implement  the  algorithm. 
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SECTION  II 

THE  EXTENDED  KALMAN  FILTER 


The  extended  Kalman  filter  is  an  extension  of  optimal 
(minimum-variance)  linear  filtering  theory  to  problems  which 
involve  significant  nonlinearities,  e.g.,  the  projectile  aero- 
dynamic coefficient  estimation  task.  This  section  provides 
the  background,  mathematical  models,  and  notational  conven- 
tions needed  for  understanding  the  EKF  algorithm  developed  in 
this  study.  The  discussion  assumes  a basic  familiarity  w'ith 
random  variables  and  state-space  notation;  additional  details 
can  be  found  in  References  1 through  8. 


2.1  KALMAN  FILTER  EQUATIONS 

To  apply  Kalman  filtering  theory  to  any  estimation 
problem,  it  is  necessary  to  derive  a linear,  stochastic,  first- 
order,  vector  matrix  differential  equation  which  models  the 
manner  in  which  the  system  states  interact  and  propagate  as  a 
function  of  time.  For  linear  systems,  this  equation  has  the 
general  form 

x(t)  = F(t)x(t)  + G(t)w(t)  + u(t)  (1) 

where  x ( t ) is  an  nxl  column  vector  representing  the  system 
state,  F(t)  is  an  n*n  dynamics  matrix  which  defines  the  inter- 
action of  the  state  vector  components,  and  v^(t)  is  a p*l 
column  vector  of  white  gaussian  noise  inputs  such  that* 

E [ w( t ) ] = 0;  Cov[w]=  E[w(t)w(i)T]  = Q(t)6(t-t)  (2) 

The  matrix  G(t)  is  an  n*p  distribution  matrix  which  indicates 
how  each  component  of  w( t ) affects  each  component  of  the  system 
state  derivative,  and  u(t)  is  a nxl  column  vector  of  known 
system  inputs.  Note  that  the  F,  G,  and  Q matrices  may  be 
time-varying.  For  a projectile  model,  the  elements  of  the 
state  vector  x will  typically  be  projectile  positions  and 
velocities,  the  elements  of  w will  be  random  inputs  such  as 
wind  disturbances  and  turbulence  effects,  and  the  elements  of 
u will  be  known  inputs  such  as  average  wind  velocity. 


♦The  symbol  E[  ] denotes  mathematical  expectation;  Cov  [w] 
denotes  the  covariance  matrix  of  w. 
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At  discrete  instants  of  time,  tR,  it  is  assumed  that 

measurements  of  linear  combinations  of  the  state  variables 
are  made.  The  equation  describing  this  measurement  process 
has  the  general  form 


*k  = Hk*k  + ^k 


(3) 


where  zR  is  a vector  of  r measured  quantities  at  time  tR,  11^ 

is  an  r*n  observation  matrix  describing  the  linear  combina- 
tions of  state  variables  which  comprise  zR  in  the  absence  of 

noise,  and  vR  is  an  r vector  of  zero  mean  gaussian  measurement 

errors  with  a covariance  matrix,  Rk,  defined  by 


(4) 


At  any  time  t,  the  objective  of  optimal  estimation 
theory  is  to  process  all  the  measurements  taken  up  to  that 
time  and  produce  an  estimate  £ ( t ) of  the  system  state  x ( t ) 
having  minimum  error,  in  a statistical  sense.  The  optimiza- 
tion criterion  most  often  chosen  is  that  of  minimizing  the 
mean  square  estimation  error.  This  minimum  mean-square  error 
estimate  is  calculated  with  the  Kalman  filtering  algorithm. 


As  measurements  become  available,  there  is  essentially 
an  instantaneous  change  in  the  knowledge  of  the  state  x(tR). 

Denoting  the  optimum  estimate  of  x(tR)  just  prior  to  the 

availability  of  zR  as  £R(-)  and  the  optimum  estimate  of  the 

state  vector  immediately  after  processing  zR  as  £R(+),  the 

Kalman  filter  generates  the  optimum  estimate  of  the  system 
state  according  to  the  following  algorithm:* 


x(  t ) = 

F ( t )x( t ) + u(t)  ; x(tR_1)  = xR_1(+) 

« % r - 1 

(5) 

*k<*> 

xk(-)  + Kk^zk  - HRxk(-)J 

(6) 

♦Only  the  continuous  form  of  the  Kalman  filter  with  discrete 
measurements  is  considered  here . 


•1 


whore  Equation  (5)  is  used  to  calculate  the  estimate  between 

measurements  and  Equation  (6)  is  used  to  update  the  estimate 

when  new  data  is  received  at  each  time  t,  . 

k 

The  n*r  matrix  Kk  is  the  Kalman  gain  matrix.  Let 
x( t ) denote  the  error  made  in  estimating  x(t),  i.e., 


x(t)  = x(t)  - x(t ) 


and  let 


P(t)  = Cov  x(t)  = E x(t)x(t)T 


Kk  is  then  computed  using  the  following  equations: 
P(t)  = F(t)P(t)  + P(t)F(t)T  + G(t)Q(t)G(t)T 
with  P(tk_1)  = Pk_1(+)  and 


Kk  ■ PkO)IlJ[HkPk<-)Hj;  ♦ fy 
Pk<  + > - [l  - XkHk]pk(-) 


where  Pk(~)  is  P(t)  just  before  the  measurement  at  time  tk  an 
Pk(+)  is  P(t)  just  after  tk>  Equation  (9)  is  used  to  calcu- 
late1 the  estimation  error  covariance  between  the  measurements 
Equations  (10)  and  (11)  are  used  to  calculate  the  Kalman  gain 
matrix  for  use  in  Equation  (6)  and  to  update  the  estimation 
error  covariance  matrix  when  a measurement  arrives. 

Figure  1 illustrates  the  structure  of  the  optimal 
linear  Kalman  filter.  This  estimation  algorithm  has  two  dis- 
tinct phases.  Equations  (5)  and  (9)  describe  the  time  evolu- 
tion of  the  state  estimate  and  its  error  statistics  between 
measurements  under  the  influence  of  system  dynamics  and  noise 
This  process  is  commonly  referred  to  as  extrapolation. 
Equations  (6),  (10),  and  (11)  indicate  how  the  estimate  and 
its  error  covariance  are  updated  at  the1  measurement  time  to 
reflect  the  new  information  available.  The  algorithm  is  opti 
mum  in  the  minimum  mean-square  error  sense  as  long  as  the 
assumed  mathematical  model  for  the  system  is  accurate. 


INITIAL  ESTIMATION 
ERROR  COVARIANCt 

P(t0) 


STATE  VARIABLE 
MODEL  OF  SYSTEM 


MEASUREMENT 

DATA 


INITIAL  STATE 
VARIABLE  ESTIMATE 

Figure  1.  Structure  of  an  Optimal  Linear  Kalman  Filter 


A unique  feature  of  the  Kalman  filter  is  that  the 
performance  analysis  is  inherent  in  the  algorithm  for  K,  . 

K 

The  matrix  P(t)  is  a complete  description  of  the  second- 
order  error  statistics  of  the  estimate.  In  particular,  the 
diagonal  terms  of  P(t)  represent  the  minimum  mean-square 
error  obtained  in  estimating  each  component  of  x(t).  Note 
P(t)  is  specified  for  all  time  by  Equations  (9),  (10),  and 
(11).  Knowledge  of  neither  x ( t ) , x(t),  nor  is  required 

to  obtain  a performance  analysis  for  the  optimal  filter. 

In  other  words,  the  performance  of  the  filter  is  completely 
determined  by  its  mathematical  model,  assuming  that  this 
model  accurately  represents  the  system  which  generates  the 
measurement  data. 

In  summary,  the  following  conditions  must  be  met  to 
implement  an  optimum  Kalman  filter: 

• The  system  must  be  linear. 

• The  matrices  F(t),  G(t),  and  u(t) 
must  be  known  functions  of  time. 

• The  vector  input  w( t ) must  be  a zero 
mean  gaussian  white  noise  process 
with  known  covariance  matrix, 

Q( t )6(t-t ) . 

• The  measurements  must  obey  Equation 
(3),  and  must  be  known  for  all  k. 


• Th«  measurement  errors  must  be  a 

gaussian  white  sequence  with  its  co- 
variance  matrix,  R.  , and  its  mean 
known . 

• To  initialize  the  filter  equations, 
the  initial  statistics  of  x must  be 
known . 

If  the  aerodynamic  coefficient  estimation  problem 
could  be  put  into  a form  which  met  all  of  the  conditions 
listed  above,  the  design  of  an  optimal  estimator  would  involve 
only  the  direct  implementation  of  the  Kalman  filter  equations. 
However,  the  projectile  dynamics  considered  here  are  nonlinear; 
the  linearity  requirement  is  violated.  Furthermore,  the  ob- 
jective of  estimating  the  projectile  aerodynamic  coej  f ic  ients , 
which  is  tantamount  to  estimating  parameters  of  the  matrix  F 
in  Equation  (1),  introduces  additional  nonlinearity.  One 
means  of  overcoming  these  problems  is  the  extended  Kalman 
filter  described  in  the  next  subsection. 


2.2  EXTENDED  KALMAN  FILTER  EQUATIONS 

Since  the  problem  under  consideration  cannot  be 
realistically  modeled  as  a linear  system,  a nonlinear  estima- 
tion technique  must  be  used.  One  method  is  the  extended 
Kalman  filter  which  is  essentially  a conventional  Kalman 
filter  design  applied  to  a mathematical  model  of  the  system 
obtained  by  linearizing  the  system  about  the  current  state 
estimate.  The  structure  of  this  algorithm  is  illustrated  in 
Figure  2.  Note  that,  because  of  the  linearization  procedure, 
the  covariance  calculation  is  now  dependent  upon  the  state 
estimate.  Consequently,  it  is  not  possible  to  calculate  the 
covariance  matrix,  as  a function  of  time,  off-line  since  it 
is  dependent  upon  the  measurement  data.  The  extended  Kalman 
filter  yields  very  nearly  optimal  estimates  if  the  lineariza- 
tion is  accurate,  i.e.,  as  long  as  the  state  estimate  is  close 
to  the  true  system  state. 

A reasonably  general  mathematical  rrodel  for  nonlinear 
stochastic  systems  is  given  by  the  equations 

x(t)  = ltx(t),  t)  + G(t)w(t)  (12) 

*k  - Hk[x(tk)]  + ; k = 1,2,...  (13) 
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Figure  2.  Structure  of  an  Extended  Kalman  Filter 


where  jf  and  are  nonlinear  differentiable  functions  of  the 
state  vector  x,  and  w(t)  and  are  zero  mean,  independent 

gaussian  white  noise  processes  having  spectral  density  and 
covariance  matrices,  Q(t)  and  R^,  respectively.  The  measure- 
ments z k are  taken  at  discrete  times  t^. 

The  first  approach  one  might  use  to  derive  a 
filtering  algorithm  for  x ( t ) in  Equation  (12)  is  to  linearize 
the  nonlinear  functions  f_  and  h^  about  an  appropriate  known 

reference  trajectory  x(t),  and  then  apply  conventional  linear 
estimation  theory,  i.e.,  the  Kalman  filter  discussed  in  the 
last  subsection.  Thus,  denoting  x(tk)  by  xk'  the  expressions 


f(x,  t)  = f ( x , t)  + 


3f 


9x 


(x-x) 


x=x 


(14) 


3h, 


bk^ik)  * ilk^k5  + 


(15) 


— k~— 


may  be  substituted  into 
corresponding  Kalman  fi 
x , 6x  ( t ) = x( t ) - x(t), 
the  reference  trajector 
estimate  of  the  state, 
as  an  extended  Kalman  f 
the  latter  are  given  in 
matrix  P(t)  is  a first- 


Equations  (12)  and  (13)  to  derive  the 
Iter  which  estimates  the  variation  in 
from  the  reference  trajectory.  When 
y is  chosen  to  be  the  current  best 
x ( t ) , the  resulting  algorithm  is  known 
ilter;  the  mechanization  equations  for 
Table  1 (see  also  Reference  1).  The 
order  approximation  to  the  estimation 
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TABLE  1. 


SUMMARY  OF  THE  EXTENDED  KALMAN  FILTER 


System  Model 
Measurement  Model 

Initial  Conditions 
Other  Assumptions 

State  Estimate  Propagation 
Error  Covariance  Propagation 

State  Estimate  Update 
Error  Covariance  Update 


error  covariance  matrix,  and  xk(-)  and  Pk(-)  denote  the  solu- 
tions to  the  propagation  equations  at  time  t.  just  before  the 

t h ^ 

k measurement  is  processed.  The  principal  practical  dif- 
ference in  mechanizing  the  extended  and  conventional  Kalman 
filters  is  that  the  gains  K^  for  the  former  depend  upon  the 

estimate;  therefore,  K^  must  be  computed  online.  Consequently, 

the  computational  burden  of  the  extended  filter  is  greater. 

Note  that  the  equations  in  Table  1 reduce  to  the  optimal 
Kalman  filter  outlined  in  the  last  section  if 

f(x(t),  t)  = F(t)x(t)  + u ( t ) (16) 

and 

= ,Ik-(tk)  (17) 
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As  a practical  matter,  one  of  the  most  important 
aspects  of  the  EKF  is  the  accuracy  of  the  linearization 
(Equations  (14)  and  (15))  about  the  state  estimate  £(t).  If 
the  estimation  error  is  large,  this  linearization  is  poor  and 
the  filter  may  not  operate  correctly.  However,  for  the  appli- 
cation considered  here  it  has  been  found  that  the  assumed 
linearization  is  satisfactory  if  the  EKF  is  properly  initial- 
ized. The  details  of  the  initialization  procedure  are  given 
in  Section  III  where  the  EKF  equations  are  applied  to  the 
aerodynamic  coefficient  estimation  problem. 
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SECTION  III 


APPLICATION  OF  THE  EXTENDED  KALMAN  FILTER 


This  section  describes  how  the  EKF  equations  in  sub- 
section 2.2  are  applied  to  aerodynamic  coefficient  estimation 
The  equations  of  motion  of  the  system  model  are  given,  and  an 
EKF  design  is  developed.  Performance  results  are  given  in 
Section  IV. 


t 


3.1  EQUATIONS  OF  MOTION  - THE  SYSTEM  MODEL 

This  subsection  presents  the  system  model  upon  which 
the  extended  Kalman  filter  design  used  in  this  study  is  based. 
The  problem  under  consideration  is  that  of  a six-degree-o f- 
freedom  rotationally  symmetric  rigid  body  flying  through  a 
ballistic  test  range.  The  projectile  dynamic  model  is  taken 
directly  from  Reference  7.  The  equations  of  motion  are  de- 
rived in  a fixed-plane  coordinate  system  and  assume  that  the 
aerodynamic  coefficients  are  expanded  as  polynomial  functions 
of  the  sine  of  the  total  angl e-of-attack . 

Figure  3 illustrates  the  two  coordinate  systems  of 
interest.  The  first  (x,y,z)  is  a fixed  inertial  system  which 
assumes  a flat,  nonrotating  earth.  The  second  (x',y',z')  is 
a fixed-plane  axis  system  attached  to  the  projectile  but 
having  zero  roll  angle.  In  this  system,  the  x'  axis  lies 
along  the  projectile  axis  of  symmetry,  and  the  origin  of  the 
system  is  fixed  to  the  projectile  center-of-gravity . The 
(x',y',z')  system  is  obtained  by  rotation  of  the  (x,y,z)  sys- 
tem through  the  two  Euler  angles  ip  and  0 in  the  indicated 
sequence.  The  projectile  roll  angle  is  measured  clockwise 
looking  downrange,  i.e.,  from  the  tail  of  the  projectile. 
Figure  4 illustrates  the  Euler  angle  rotations  and  the  rela- 
tionship between  the  two  coordinate  systems. 

Twelve  state  variables  are  used  to  define  the  six- 
degree-of-f reedom  dynamic  equations  of  the  projectile.  These 
state  variables  are  summarized  in  Table  2 with  their  defini- 
tions and  units.  Table  3 defines  the  constants  required  by 
the  model,  Table  4 defines  the  intermediate  variables  used, 
and  Table  5 gives  the  polynomial  expansions  of  the  aerodynamic 
coefficients  which  are  to  be  estimated. 
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TABLE  2.  SYSTFM  STATE  VARIABLES 


SYMBOL 

DEFINITION 

UNITS 

X 

Downrange  position 

1 ft 

y 

Crossrange  position 

z 

Vertical  position 

ft 

u 

Velocity  along  x’  axis 

! ft /sec 

v 

Velocity  along  y’  axis 

f t / sec 

w 

Velocity  along  z’  axis 

f t/sec 

■p 

Euler  angle) 

fSee  Figure  4 

rad 

0 

Euler  angle’ 

rad 

■t> 

Roll  angle  about  x axis 

rad 

Time  derivative  of  ’p 

rad/sec 

• 

9 

Time  derivative  of  9 

rad/sec 

P 

Spin  rate  about  x'  axis 

rad/sec 

TABLE  3.  SYSTEM  CONSTANTS 


i 


t 


SYMBOL 


DEFINITION 


UNITS 


ft 


ft' 


Reference  diameter 

2 

Reference  area;  A = — — 

4 

Axial  moment  of  inertia 
Transverse  moment  of  inertia 
Mass 

Reference  velocity 

Acceleration  of  gravity,  g=32. 17405 
if  = 3.1415926536 


d 

A 


y 

m 

/ 

o 

g 


slug-f t' 
* 

slug-f t' 
slugs 
f t/sec 
ft /sec 


* 


» TABLE  4.  SYSTEM  INTERMEDIATE  VARIABLES 


SYMBOL 

DFFINITION 

UNITS 

Air  ci  it}-; 

„r  T 4. 2561 

. = 2.  _<7692x’.  0--3|l-(6.S754.\10-6)z 

lb-s;<  2 ' ft4 

v 

Velocity  ‘luinttude  : 

\ = \ 11  - + 

f t ■ sec 

Dynamic  pressure; 

■1 

q - 1/2  V2 

9 

lb  ft“ 

Sine  of  the  total  angle-of-attack ; 
e **  \v*“  ♦ w2/v 

TABLE  5.  AERODYNAMIC  COEFFICIENTS  AND  THEIR 
POLYNOMIAL  EXPANSIONS 


* 

SYMBOL 
■ — ■ ■ 

DEFINITION  AND  EXPANSION* 

If 

c,. 

Axial  force  coefficient, 

r 

S'  * S + CN2  2 * CXV  <VV) 

Sci 

Normal  force  coefficient  slope; 

i 

S,  - CN1  + CN,3-2 

> 

CVpa 

Magnus  force  coefficient  slope, 

^Ypi  = CYp  a + CYpa3c2 

C 

mi 

Pitching  moment  coefficient  slope 
i Cma  + ^mi3"  + ^maV 

i 

• 

C 

mq 

Damping  moment  coefficient; 

C - C ♦ c „c2 
mq  mq  mq2 

C 

npi 

Magnus  moment  coefficient  slope 

c « c * c • 2 

np t npi  npi3~ 

ciP 

Spin  deceleration  coefficient 

<v 

Fin  cant  moment  f \c  lent 

ci 

Moment,  of  in«  r*m  ratio  coefficient;  C.  * 1 

*c  is  the 

sine  oi  the  total  angle-of -at tack . 

i 


Given  the  variables  and  constants  defined  in  Tables 
2 through  5,  the  twelve  state  dynamic  equations  can  be  written 
as  shown  in  Table  6.  The  moment  of  inertia  ratio  coefficient, 
Cj , used  in  Equation  11  in  Table  6 is  equal  to  one  in  any 

real  system.  However , the  filter  can  estimate  the  ratio  I /I 

x'  y 

more  accurately  than  it  can  be  measured.  Since  a measured 
value  of  Ix/I  is  used  in  the  filter  implementation,  the  esti- 
mate of  Cj  adjusts  for  the  error  in  this  measured  quantity. 


TABLE  6.  EQUATIONS  OF  MOTION 


EQUATION 


= u co.s0  cosiji  - v simli  + w sin0  rosty 
' = u cos 9 si nii  + v cosi|i  + w sin'1  sinv 


. = -u  sinfl  + w cos't 


u = C^.  + r sin1-  - u + ,v  cos') 

1 ’ -ft?)  ? * ($)($)  V.  v - «.  • 

* - -ft?)  s,  ? - ft?)(i)  'V,,,  s - « * 


( ii  cos  * ■ + w sin9) 


v sin1  * n ii 


5 = |-(  $ ♦ (^)(&)  % • * f^)®)  cnpt  ? 

* ^ sinn  + ri  (!")}?.7,V- 

5 = (t5)  v v * (¥)(&)  % • * (¥-)(£-?)  (H|1'«  ; 

' y «.*  • y / » v f 


P cos)  Cjl  -- J - . cos-  sin* 


• (^-)(£)  <V  • (t‘)  <■ 
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This  completes  the  definition  of  the  system  dynamic 
model.  The  model  is  used  to  generate  projectile  trajectories 
from  which  measurements  are  derived  to  be  applied  to  the 
filter;  it  is  also  used  as  the  mathematical  model  upon  which 
the  filter  design  is  based. 

The  measurement  model  takes  noisy  measurements  of 
position  and  angular  orientation  (x,y,z,^,e,  and  <t> ) and  time 
as  the  projectile  passes  50  downrange  measurement  stations 
during  its  flight.  To  each  of  these  quantities,  a random 
gaussian  uncorrelated  noise  sequence  is  added  with  a given 
rms  value. 


3.2  EXTENDED  KAI.MAN  FILTER  ALGORITHM  DESIGN 


In  this  section,  the  EKF  algorithm  summarized  in 
subsection  2.2  is  applied  to  estimate  the  17  aerodynamic 
coefficients  contained  in  the  system  model  defined  in  sub- 
section 3.1.  To  do  this,  the  system  model  must  be  put  into 
the  form  of  Equations  (12)  and  (13).  Let  the  12  dynamic  state 
variables  listed  in  Table  2 be  brought  together,  in  the  order 
listed,  into  the  vector  and  let  the  17  constant  aerodynamic 
coefficients  form  the  vector  a..  Then  the  equations  of  motion 
listed  in  Table  G can  be  expressed  in  state  vector  notation 
as 


(18) 


where  £ is  a nonlinear  vector  function  of  £ and  a.  Since  the 
unknown  coefficients  contained  in  a_  are  assumed  to  be  con- 
stant, the  dynamic  model  for  a_  is  given  by 


a = 0 


(19) 


where  0 is  the  zero  vector. 


To  apply  the  EKF  algorithm,  all  of  the  quantities  to 
be  estimated  in  Equation  (18)  must  be  expressed  in  state 
variable  form.  This  is  readily  accomplished  by  augmenting 
Equation  (18)  with  Equation  (19)  as  follows: 


1 

~£(£.a)" 

a 

0 

Equation  (20) 


is  a special  case  of  Equation  (12)  where 


x = 

V 

; f ( x , 1 ) * 

‘g(Z-a)" 

a 

0 

— — 
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G(t)w(t)  = 0 ( 22  j 

The  measured  quantities  in  this  application  are  sys- 
tem states  except  for  the  time  measurement.  Thus,  the  state 
measurements  are  linear,  and  Equation  (13)  can  be  expressed 


= + v. 


k=l,2 50 


where  tk  is  the  measured  measurement  time  and  z k are  the  noisy 

measured  values  of  x,  y,  z,  i p,  9,  and  <t> . For  this  case,  the 
matrix  is  time-invariant  and  given  by 

fl  00000000 


6*20  j (24) 


The  vector  Vj,  in  Equation  (23)  represents  the  measure- 
ment noise  and  consists  of  the  combined  effect  of  the  position 
and  time  measurement  errors.  Since  the  random  time  measure- 
ment error,  c1 , is  very  small  for  this  application  (0.5  ysec 

rms),  the  projectile  velocities  can  be  assumed  constant  over 
the  error  interval  and  the  total  measurement  error  is  approxi- 
ma  t ed  by 

vk  = v ’ - Hkx(tk)Et  (25) 

The  covariance  of  vk  is  thus  given  by 


Rk  ' Rk  + Hk  E[i<vi(tk>T]Hj  °2t 


where  v^  is  the  measurement  error  without  time  measurement 

error,  its  covariance  matrix,  and  the  rms  value  of  the 

time  measurement  error.  For  the  current  application,  the 
second  term  on  the  right-hand  side  of  Equation  (26)  is  sig- 
nificant only  in  the  measurement  of  the  projectile  roll  angle, 
which  changes  at  a high  angular  rate.  Consequently,  Rk  is 

assumed  equal  to  Rk  except  for  the  diagonal  element  associated 
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with  the  roll -angle  measurement  error,  which  is  approximated 
by 

Rk(6’6)  £ Rk  (6’6)  + P2  °?  (27) 

where  p is  the  estimate  of  projectile  spin  rate.  The  validity 
of  this  approximation  is  demonstrated  by  the  results  presented 
in  Section  IV  of  this  report. 


Before  presenting  the  performance  and  sensitivity 
results,  a few  specific  implementation  details  should  be  dis- 
cussed. First,  the  filter  design  is  based  upon  the  correct 
values  of  the  projectile  constants,  d,  Ix,  I m,  and  g. 

Normally  each  of  these  quantities  will  be  measured  off-] ine 
with  some  errors.  Except  for  the  ratio  I /I  , none  of  the 

t ^ y 

projectile  constants  can  be  estimated  by  the  filter  indepen- 
dently from  the  aerodynamic  coefficients.  For  example,  in 
Equation  4 of  Table  6,  is  estimated  by  assuming  that  the 

value  of  A/m  is  known.  In  practice,  the  filter  can  estimate 
only  the  combination,  A/m  ?x>  not  Cx  and  A/m  individually. 

In  this  sense,  the  problem  is  over  parameterized,  and  all 
measured  or  unknown  quantities  cannot  be  estimated  indepen- 
dently. Therefore,  for  convenience  any  measurement  errors  in 
the  projectile  parameters  are  attributed  to  uncertainty  in 
, the  aerodynamic  coefficients,  in  order  to  maintain  the  expli- 

cit form  of  the  equations  of  motion  in  Table  6. 

Second,  the  assumption  of  a flat  non-rotating  earth 
is  not  valid  for  data  generated  in  a ballistic  range  and  must 
be  accounted  for  in  an  operational  program.  Third,  the  compu- 
tation of  the  F(x(t),t)  matrix  in  Table  1 requires  the  compu- 
' tation  of 

9 _f  ( x , t ) 

. F(x,t)  = (28) 

9x(  t ) 

This  matrix  of  partial  derivatives  is  explicitly  derived  for 
the  filter's  implementation.  The  only  assumption  made  in 
this  derivation  is  that  the  air  density,  p,  is  independent  of 
altitude,  z.  While  not  strictly  correct,  the  error  induced 
by  this  approximation  is  trivial.  Appendix  A gives  the  de- 
tailed equations  used  to  compute  this  matrix. 


Finally,  the  initialization  procedure  is  very  impor- 
tant to  the  proper  operation  of  the  filter.  Before  the  first 
, measurement,  only  very  rough  a priori  estimates  of  the 
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projectile's  position  and  downrange  velocity  at  the  first  mea- 
surement station  are  available.  Thus,  the  angular  and  trans- 
lational position  and  velocity  initial  estimates  are  set  to 
zero  or  nominal  values  for  this  station,  and  their  error  co- 
variances are  initialized  at  correspondingly  high  values. 

The  coefficient  estimates  and  their  error  covariances  are 
initialized  based  upon  the  a priori  knowledge  of  the  projec- 
tile's aerodynamic  properties.  Based  on  these  initial  con- 
ditions, the  first  measurement  set  is  processed,  and  the  filter 
estimate  and  its  covariance  are  updated.  This  single  measure- 
ment of  six  positions  brings  the  error  covariance  of  the  posi- 
tion estimate  down  to  the  measurement  noise  level  but  provides 
no  information  on  the  six  velocities  at  the  first  station. 

Thus,  the  velocity  estimates  remain  extremely  poor. 

The  filter  estimates  and  covariance  are  propagated 
to  the  next  measurement  time,  and  the  second  filter  update  is 
performed  using  the  second  measurement  set.  This  greatly 
improves  the  estimates  of  the  six  projectile  velocities  be- 
cause two  sets  of  accurate  position  measurements  at  different 
times  have  been  processed.  However,  the  error  covariance  at 
this  point  has  not  been  computed  correctly.  Recall  that  the 
EKF  design  requires  a linearization  about  the  state  estimate 
to  propagate  the  error  covariance  between  two  measurements. 
Between  the  first  and  second  measurements  the  velocity  esti- 
mates are  very  inaccurate,  resulting  in  a poor  linearization. 
Consequently,  the  covariance  after  the  second  measurement  is 
not  a good  reflection  of  the  true  error  covariance.  This 
problem  must  be  corrected  at  the  second  measurement  if  future 
measurements  are  to  be  processed  correctly. 

The  following  filter  resetting  procedure  was  found 
to  be  a satisfactory  method  of  correcting  the  problem  de- 
scribed above.  Immediately  after  the  second  measurement  up- 
date, the  filtering  problem  is  essentially  restarted  by  re- 
setting the  error  covariance  matrix  as  follows:  the  error 

covariance  matrix  off  diagonal  elements  are  all  set  to  zero; 
the  diagonal  elements  corresponding  to  positions  are  reset  to 
the  measurement  noise  levels  for  those  positions;  the  diagonal 
elements  corresponding  to  coefficient  estimates  are  set  to 
their  original  values  assumed  before  the  first  measurement; 
the  diagonal  element  corresponding  to  the  velocity  along  the 
projectile  axis  of  symmetry  is  unchanged;  and  the  diagonal 
elements  corresponding  to  the  other  two  translational  veloci- 
ties and  the  three  angular  velocities  are  increased  by  a fac- 
tor of  ten  over  their  current  values.  This  covariance  matrix 
is  then  propagated  together  with  the  estimate  to  the  third 
measurement,  the  third  update  occurs,  and  normal  processing 
continues  thereafter.  Since  the  estimates  of  position  and 
velocity  are  fairly  accurate  after  the  reset,  the  lineariza- 
tion is  good  over  the  propagation  between  the  reset  and  the 


third  measurement  and  the  filter  operates  properly,  correc-tl 
computing  the  estimates  and  their  corresponding  error  covari 
ance  matrix.  Essentially,  the  reset  starts  the  filter  over 
after  the  second  measurement  update  with  a new  initial  co- 
variance  matrix  that  is  a much  more  accurate  representation 
of  true  velocity  estimation  errors.  The  effectiveness  of 
this  initialization  procedure  is  demonstrated  by  the  experi- 
mental results  presented  in  Section  IV. 
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SECTION  IV 


PERFORMANCE  AND  SENSITIVITY  RESULTS 


This  section  presents  performance  and  sensitivity 
results  obtained  using  the  EKF  algorithm  described  in  sub- 
section 3.2.  Test  procedures  and  performance  measures  are 
discussed,  a nominal  test  case  is  defined,  and  the  filter 
performance  for  this  case  is  demonstrated.  Finally,  sensi- 
tivities to  variations  in  the  nominal  test  conditions  are 
invest igated . 


4.1 


TEST  PROCEDURES  AND  PERFORMANCE  MEASURES 


Experimental  data  is  generated  by  a trajectory  simu- 
lation program  which  integrates  the  equations  of  motion  given 
in  Table  6,  based  on  given  true  values  for  all  the-  components 
of  the  model.  As  the  downrange  position  of  the  projectile 
reaches  each  measurement  station,  the  time,  position,  and 
angular  orientation  of  the  projectile  are  computed  and  cor- 
rupted by  independent  zero-mean  gaussian  noise  sequences. 

Thus,  the  input  data  to  the  filter  consists  of  50  sets  of 
seven  measurements  (three  angular  positions,  three  transla- 
tional positions,  and  time).  The  filter  design  is  specified 
by  the  system  model  and  its  constants.  It  begins  with  a given 
initial  state  estimate  x„  and  an  assumed  covariance  matrix, 

P ~ 1 u 

1 0 • 

The  performance  of  the  filter  is  evaluated  by  observ- 
ing its  estimation  error  in  relation  to  the  error  standard 
deviations  obtained  from  the  covariance  matrix  computation. 

If  the  errors  are  consistent  with  the  computed  standard  devia- 
tion, then  the  filter  is  considered  to  be  operating  properly 
and  its  absolute  performance  can  be  judged  by  the  estimation 
errors  achieved.  Thus,  there  are  two  basic  measures  of  per- 
formance, estimation  error  and  consistency  between  this  error 
and  its  computed  standard  deviation. 


4 . 2 


THE  NOMINAL  CASE 


The  nominal  trajectory  is  obtained  by  integration  o f 
the  equations  of  motion  in  Table  G,  using  the  projfjctile  con- 
stants listed  in  Table  7,  the  trajectory  dynamic  initial  con- 
ditions (at  t=0)  listed  in  Table  8,  and  the  true  aerodynamic 
coefficient  values  listed  in  Table  9.  As  the  downrange  posi- 
tion (x)  of  the  projectile  passes  the  50  measurement  station 
locations  listed  in  Table  10,  measurements  of  time,  angular 
position,  and  translational  position  are  made.  These 
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TABLE  7. 

NOMINAL  PROJECTILE  CONSTANTS 

SYMBOL 

DESCRIPTION 

NOMINAL  VALUE 

d 

Diameter 

9 . 8333xl0-2  ft 
( 3 3 mm ) j 

m 

Mass 

2.4865xl0-2  slug 
(0.80  ff)* 

I 

X 

Axial  moment  of  inertia 

3.2376xl0~5  slug-ft2 
(0.15  tt  in  . 2 ) 

:y 

Transverse  moment  of  inertia 

2.6764xl0"4  slug-ft2 
(1.24  P in. 2) 

vo 

Reference  velocity 



3 . 34 57xl03  ft/sec 
(Mach  3.00) 

*32.174  ft  = 1 slug 


TABLE  8.  TRAJECTORY  GENERATOR  AND  FILTER 
DYNAMIC  STATE  INITIAL  CONDITIONS 


ft 

0 

.0 

1 .31  x 

10-4 

0 

1 

.0 

n 

20 

0 

20 

.008 

20 

.0 

1 

,11 

ft/  isei: 

:t , :t5 1 

,0 

3,34  0 

.7 

3,000 

. t) 

600 

.0 

f 

" 

.0 

-201 

.0 

0 

.0 

300 

.0 

ft/  M;i' 

n 

.0 

1 37 

.2 

0 

.0 

300 

.0 

rad 

0 

. 1) 

0 

. 000 

0 

.0 

1 

.0 

rad 

-1  .745 

x 10-'* 

0 

.04  0 

0 

.0 

1 

,0 

rad 

0 

, 0 

10 

.05 

0 

.0 

200 

.0 

rad/soc 

0 

,0 

50 

. K5 

0 

.0 

100 

. () 

rad/sf< 

55 

.0 

-IK 

. 1 1 

0 

. 0 

100 

.0 

rad/sia: 

10,703 

,u 

10,7oo 

. K 

10,000 

.0 

1 ,000 

.0 

If 


AERODYNAMIC 

COEFFICIENT 

TRUE  VALUE 
— 

INITIAL 

FILTER 

ESTIMATE 

RMS  ERROR 
IN  INITIAL 
ESTIMATE 
ASSUMED  BY 
FILTER 

INITIAL 

PERCENT 

ERROR 

c* 

0.225 

0.300 

0.2 

33 

c 

\X2 

1.5 

0.0 

1.0 

-100 

r 

^XV 

-0.54xl0-4 

0.0 

1 .0*10~4 

-100 

Cv 

Na 

2.87 

2,75 

0.32 

-4 

CNa3 

10.0 

0.0 

5 . 0 

-100 

CYpa 

-0.9 

-1.0 

0.1 

11 

Spa  3 

5.0 

0.0 

3.0 

-100 

Cma 

3.15 

2.50 

0.75 

-21 

C „ 
ma3 

-6.0 

0.0 

6.0 

-100 

maV 

2.58xl0-4 

0.0 

3 . Ox 10~4 

-100 

C 

mq 

-18.0 

-15.0 

5.0 

-17 

Crr.q2 

10.0 

0.0 

5.0 

-100 

C 

npa 

0.3 

0.0 

1.0 

-100 

c 

npa3 

2.0 

0.0 

2.0 

-100 

% 

-0.024 

-0.015 

0.015 

-38 

cu 

-l.OxlO"3 

0.0 

_3 

1.0x10  J 

-100 

_i_i 

1.00 

1.01 



0.01 

1 
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TABLE  10. 


MEASUREMENT  STATION  LOCATIONS 


I 


STATION 

NUMBER 

DOWNRANGE  POS I T ION 
(FT) 

STATION 

NUMBER 

DOWNRANGE  POSITION 
(FT) 

1 

5 

26 

275 

2 

15 

27 

290 

3 

25 

28 

305 

4 

35 

i 

29 

320 

5 

45 

30 

335 

6 

55 

31 

350 

7 

65 

32 

380 

8 

75 

33 

395 

9 

85 

34 

410 

10 

95 

35 

425 

11 

115 

36 

440 

12 

125 

37 

445 

13 

135 

38 

470 

14 

145 

39 

485 

15 

155 

40 

500 

16 

170 

41 

515 

17 

180 

42 

530 

18 

190 

43 

545 

19 

210 

44 

560 

20 

220 

45 

585 

21 

230 

46 

600 

22 

240 

47 

615 

23 

250 

48 

630 

24 

260 

49 

645 

25 

270 

50 

660 
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measurements  consist  of  the  true  values  of  these  quantities 
corrupted  by  independent  sequences  of  uncorrelated  zero-mean 
gaussian  noise  samples.  The  standard  deviations  of  the  noise 
sequences  for  the  measurement  of  position  (translational  and 
angular)  and  time  are  given  in  Table  11.  Figures  5 and  6 are 
plots  of  some  of  the  states  of  the  nomiral  trajectory  to 
illustrate  its  general  form  and  frequency  content.  Figure  6 
also  includes  a plot  of  the  total  angle-of-attack,  a,  defined 
by 


a 


sin 


■(e) 


sin 


(29) 


This  trajectory  represents  a model  of  a 30mm  projectile  travel- 
ing at  approximately  Mach  3 with  initial  Euler  angle  rates  of 
Q(0)  = 55  radians  per  second  and  i|*(0)  = 0 radians  per  second. 

The  filter  design  is  based  upon  the  same  equations 
of  motion  and  projectile  constants  as  the  nominal  trajectory. 
The  filter  is  initialized  at  time  tj  when  the  projectile 

activates  the  first  measurement  station.  Table  ft  gives  the 
true  values  of  the  projectile  dynamic  states  at  t^ , the  ini- 
tial filter  estimates,  and  the  initial  rms  errors  in  the  esti- 
mates assumed  by  the  filter  at  time  t1 . Note  that  only  rough 

knowledge  of  the  initial  positions  and  velocities  is  assumed 
by  the  filter. 


TABLE  11.  NOMINAL  MEASUREMENT  ERROR  STANDARD  DEVIATIONS 


SYMBOL 

DESCRIPTION 

_ ...  . . 

NOMINAL  VALUE 

c 

— -J 

Translational  Position  (x,y,z) 

0.01  ft 

P 

Measurement  Error  Standard 
Deviation 

(-0.1  in.  ) 

0 

a 

Angular  Orientation  (^,0,*) 
Measurement  Error  Standard 
Deviation 

1 . 73x 10~3  rad 
(0.1  deg) 

ct 

Time  of  Measurement  Error 
Standard  Deviation 

0.5  usee 
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In  a similar  manner,  Table  9 gives  the  true  values 
of  the  17  aerodynamic  coefficients,  the  initial  values  of  the 
coefficient  estimates,  and  the  rms  estimation  errors  assumed 
by  the  filter.  Table  9 also  lists  the  percent  errors  in  the 
initial  coefficient  estimates.  This  completes  the  specifica- 
tion of  the  trajectory,  measurement  process,  filter  design, 
and  filter  initialization  for  the  nominal  case. 

Figures  7 through  13  illustrate  the  coefficient  esti- 
mation error  performance  of  the  extended  Kalman  filter  for 
the  nominal  case.  The  estimation  errors  (solid  lines)  are 
plotted  as  percent  errors  relative  to  the  true  coefficient 
values,  that  is, 


Percent  estimate  error  = 100 


(C  - C) 


where  C represents  the  estimate  of  a coefficient,  C its  true 
value.  Also  plotted  on  the  same  figures  (dashed  lines)  are 
the  normalized  positive  and  negative  one-o  values  of  these 
errors  as  computed  by  the  filter  covariance  matrix.  These 
one-o  values  are  defined  as 


Normal  iz.ed  computed 
standard  deviation 


= 100 


where  P^,  is  the  value  of  the  diagonal  element  in  the  filter 

covariance  matrix  associated  with  C.  Table  12  summarizes  the 
filter  errors  at  the  end  of  the  trajectory. 

Figures  7 through  13  indicate  that,  in  all  cases, 
the  parameter  estimation  errors  achieved  are  consistent  with 
the  standard  deviations  computed  by  the  filter.  Some  of  the 
coefficients  are  estimated  very  accurately  by  the  filter  while 
others  are  not  estimated  at  all.  This  is  due  to  the  fact 
that  some  of  the  coefficients  have  very  little  effect  on  the 
trajectory  and  are  therefore  not  visible  in  the  data.  Others 
strongly  affect  the  trajectory  and  are  easily  estimated.  The 
covariance  calculation  correctly  distinguishes  between  these 
coefficient  types  and  extracts  most  of  the  available  informa- 
tion from  the  data. 


An  attempt  was 
nominal  case  by  filteri 
This  was  done  using  the 
respective  covariances 
the  second  pass  through 
the  second  pass  with  be 
assumed  in  the  F.KF  deri 
the  filter  should  be  cl 
lists  the  results  of  th 


made  to  improve  the  estimates  for  the 
ng  the  nominal  data  set  a second  time. 

final  coefficient  estimates  and  their 
as  the  filter  initial  condition  for 
the  data.  Since  the  filter  begins 
tter  estimates,  the  linearization 
vation  should  be  more  accurate,  and 
oser  to  an  optimum  design.  Table  12 
is  experiment,  designated  Test  No.  1, 
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Figure  10.  Coefficient  Percent  Estimation  Errors  and  Their 

Computed  Standard  Deviations  for  C , C „ , and 
„ _ mm  mm 


Nominal  Case 


COMPUTED  STANDARD  DEVIATION 


ESTIMATION  ERROR 


MEASURE  ME  NT  NUMBER 


COMPUTED  STANDARD  DEVIATION 


ESTIMATION  ERROR 


measurement  NUMBER 


Figure  11 


Coefficient  Percent  Estimation  Errors  and  Their 
Computed  Standard  Deviations  for  C and  C 
Nominal  Case  mq  nq2 
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TABLE  112.  FILTER  PERFORMANCE  TEST  RESULTS 
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as  well  as  the  results  for  other  tests  not  yet  discussed. 

The  second  pass'"ttrrough  the  data  results  in  no  real  improve- 
ment in  estimation  accuracy,  indicating  that  the  first  pass 
extracted  almost  all  of  the  available  information. 


Another  illustration  of  the  fact  that  all  pertinent 
information  about  the  coefficients  has  been  extracted  from 
the  data  is  obtained  by  integrating  the  trajectory  using  the 
estimated  coefficients  and  comparing  it  with  the  nominal  tra- 
jectory based  on  the  true  coefficients.  Figure  14  shows  the 
trajectory  states  y,  v,  and  <l>  for  the  trajectory  generated 
using  the  estimated  values  of  the  coefficients  provided  by 
the  EKF  for  the  nominal  case.  Ttie.se  curves  are  virtually 
indistinguishable  from  their  true  counterparts  in  Figures  5 
and  6. 


0 0 Oh  0 10  O 1 r. 


States  y,  v,  and  'y  for  the  Trajectory  Based  on 
the  Estimated  Aerodynamic  Coefficients  for  the 
Nominal  Case 


PERFORMANCE 
ERROR  LEVEL 


SENSITIVITY  TO  TRAJECTORY  AND  MEASUREMENT 


Table  12.  summarizes  th*  performance  of  the  EKF  for 


the  nominal  case  and  several  additional  simulations.  Tests 
No.  2 and  No.  3 are  designed  to  show  the  sensitivity  of  EKF 
performance  to  the  amount  of  angular  motion  in  the  trajectory. 
The  ~t  ondi  t iou£_Jjr>r  these  two  tests  differ  from  the  nominal 
case  only  in  the  lnYTtsT-  -concH  t ions  on  the  trajectory  generator 
Recall  that  for  the  nominal  t ra jeft-e-r-se.  t_he  initial  Euler  angle 
rates  are  _ 


0(0)  = 55  radians  per  second 
V ( 0 ) = o radians  per  second 


In  Test  No.  2,  6(0)  is  set  to  zero  resulting  in  a trajectory 
without  the  high  frequency  nutational  motion  seen  in  Figures 
5 and  6.  In  Tost  No.  3,  j(0)  is  set  equal  to  100  radians  j < r 
second  yielding  a trajectory  with  a nutational  motion  of 
approximately  twice  the  amplitude  of  the  nominal  trajectory. 


In  each  case  the  filter  operated  in  a consistent 
manner.  The  final  estimation  errors  and  their  computed  stan- 
dard deviations  are  shown  in  Table  12.  For  Test  No.  2,  the 
lack  of  angular  motion  severely  affects  the  filter's  ability 
to  estimate  those  coefficients  associated  with  projectile 
nutations.  This  problem  is  correctly  reflected  by  the  co- 
variance  calculations.  The  lack  of  angular  motion  helps, 


in  the  estimation  of  the 


drag  coefficients  Cy 

and  in  the  spin  damping  coefficients  C, „ and  C, 


however , 

XY  mhu  til  o t * v -■  o|.-  J.  11  uuii.j/nif;  v.  wr  i i ini  l ^ -J  p ci  1 1 vi  ^ 

No.  3 has  large  amplitude  angular  motions  which  aid  in 
estimation  of  coefficients  associated  with  this  motion 


and 

Test 

t he 
For 


the 
It  is  also 


estimates  of  Cx2,  C and  C 
interesting  to  note 


ma  3 
that 


the  estimates  of  Cv  and  C, 


are  greatly  im- 
there  is  a degra- 


from  the  nominal  cast- 


example  , 
proved . 

dation  in  — ~ ^ dm* 

In  Test  No.  4,  the  time  measurement  noise  is  set  to  zero  both 
in  the  measurement  model  and  in  the  filter  computations;  other- 
wise, everything  is  the  same  as  the  nominal  case.  Time  mea- 
surement noise  apparently  has  a negligible  effect  on  the  esti- 
mation of  all  coefficients  except  those  associated  with  the 
spin  rate,  Cj  and  C^.  Due  to  the  high  spin  rate  of  about 

10,000  radians  per  second,  t ime  errors  look  like  increased 
roll  position  measurement  errors  as  described  in  subsection 
3.2.  With  no  time  measurement  noise,  the  effective  errors  in 
roll  angle  measurement  are  smaller,  thereby  reducing  the 
coefficient  estimation  errors. 


In  Test  No.  5,  the  translational  and  angular  position 
measurement  noise  variances  are  raised  by  a factor  of  ten. 

This  change  is  made  in  both  the  measurement  error  model  and 
in  the  filter  design.  As  might  be  expected,  the  higher  noise 
environment  degrades  the  parameter  estimates.  In  particular, 
the  ability  to  estimate  is  almost  lost,  and  those  coeffi- 

cients which  were  accurately  estimated  in  the  nominal  case 
are  generally  degraded  by  factors  of  two  or  more  in  the  higher 
noise  environment.  This  effect  is  to  be  expected  and  it  is 
accurately  accounted  for  by  the  FKF ' s covariance  calculation. 

In  summary,  the  study  demonstrates  that  aerodynamic 
coefficient  estimation  accuracy  is  strongly  influenced  by  the 
test  scenario  in  which  the  filter  operates.  However,  the  FKF 
algorithm  generally  extracts  most  of  the  information  about 
each  coefficient  available  from  the  data  and  correctly  assigns 
uncertainties  to  the  estimates. 
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PERFORMANCE  SENSITIVITY  TO  MODELING  ERROR 


In  subsection  4.3,  all  of  the  experiments  involve  a 
filter  design  where  the  filter  assumes  the  trajectory  model 
actually  used  to  generate  the  data.  In  many  practical  appli- 
cations, the  filter  design  model  may  only  be  an  approximation 
to  the  system  which  generates  the  data.  Such  modeling  error 
can  severely  affect  the  performance  of  the  filtering  algorithm 
Table  13  summarizes  the  results  of  the  nominal  case  and  three 
experiments  designed  to  test  the  effect  of  modeling  error  on 
filter  performance. 


TABLE  13.  FILTER  SENSITIVITY  TEST  RESULTS 
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Test  No.  6 investigates  the  influence  of  unmodeled 
time  measurement  noise.  Here  the  filter  assumes  that  the  rms 
time  measurement  errors  are  zero  when  they  are  actually  the 
same  as  for  the  nominal  case.  The  effect  of  this  modeling 
error  is  significant  only  in  the  estimates  of  those  coeffi- 
cients associated  with  the  spin  rate  dynamics,  Cj  and  C.j  v . 

While  the  table  seems  to  indicate  performance  improvement, 
this  is  not  necessarily  the  case  since  the  errors  involved 
are  random.  Figure  15  show's  the  estimation  error  and  standard 
deviation  of  C^  for  the  nominal  case.  Figure  16  shows  them 

for  Test  No.  6:  evidently  the  standard  deviation  is  lower  but 
the  actual  estimation  error  is  larger.  Thus,  the  average- 
error  performance  is  poorer  than  nominal,  and  the  standard 
deviation  gives  a less  conservative  indication  of  the  rms 
error.  The  same  effect  is  shown  in  Figures  17  and  18  for 
C-^.  Figures  19  and  20  show  the  effect  of  this  modeling  error 

on  the  filter  roll  angle  residual  process.  This  process  is 
defined  to  be  tho  difference  between  each  roll  angle  measure- 
ment and  the  value  of  that  measurement  predicted  by  the  filter 
based  upon  past  data.  In  these  figures,  the  residuals  are 
normalized  by  their  computed  standard  deviations  which  are 
calculated  from,  the  filter  covariance  matrix.  The  residual 
process  for  Test  No.  6 in  Figure  20  has  a nuch  larger  rms 
value  than  that  for  the  nominal  case  shown  in  Figure  19.  The 
fact  that  the  rms  normalized  residuals  in  Figure  20  are  so 
large  is  a good  indication  that  modeling  error  exists  in  the 
roll  dynamics  or  measurement  model. 


Figure  15.  Percent  Estimation  Error  and  Its  Computed 

Standard  Deviation;  Nominal  Case 
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Figure  20.  Normalized  Roll  Angle  Measurement  Residuals; 

Test  No.  6 - Time  Measurement  Errors  Unmodeled 

Test  No.  7 investigates  the  effect  of  not  modeling 

the  aerodynamic  coefficients  Cvv  and  C ...  which  account  for 
_ _ ay  mV 

changes  in  Cy  and  C as  velocity  changes.  Figure  21  illus- 
trates the  effect  of  this  modeling  error  on  the  estimates  of 
"X2  anr*  ^tn'(3’  n°te  the  poor  performance  and  apparent  divergence 

as  compared  to  the  nominal  case  in  Figures  7 and  10  where 
these  estimates  are  much  better  behaved.  Clearly,  these  ve- 
locity terms  must  be  modeled  by  the  fill  t»r  if  the  aerodynamic 
coefficients  change  with  velocity.  Table  13  summarizes  the 
results  of  this  experiment. 

Test  No.  8 investigates  the  effect  of  the  filter 
assuming  rms  position  and  angle  measurement  errors  three  times 
larger  than  actually  exist.  Table  13  shows  that  there  is 
some  degradation  in  performance,  but  it  is  not  as  large  a 
might  be  expected.  The  computed  standard  deviations  are 
overly  conservative  in  this  case,  which  might  be  desirable  in 
a practical  operational  program. 

1.3  DISCUSSION  OF  RESULTS 

In  the  absence  of  significant  modeling  error,  tin 
1 KF  algorithm  demonstrates  excellent  performance.  While  the 
design  is  known  to  be  suboptimal  because  of  the  non  1 1 noa r i t i es 
in  the  projectile  dynamics,  it  appears  that  most  of  t h<  sig- 
nificant information  in  the  data  is  extracted  by  the  algorithm; 
both  the  estimation  errors  and  the  measurement  residuals  are 
completely  consistent  with  the  tarns  values  predicted  by  the 
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Figure  21.  C^2  iint^  ^nri3  Percent  Estimation  Errors  and 

Their  Computed  Standard  Deviations;  Test  No.  7 
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EKF  covariance  equations.  For  the  nominal  case,  the  linear 
drag  coefficient  Cj.  is  estimated  to  within  2 percent  and  the 

drag  variation  with  velocity  C to  within  45  percent.  The 

other  force  coefficients  have  only  a minor  influence  on  the 
trajectory  and  are  not  accurately  estimated.  The  moment  co- 
efficient C is  estimated  to  within  1 percent;  C is  esti- 
mci  r mq 

mated  to  within  3 percent,  and  to  within  12  percent. 

The  spin  coefficient  Cj  is  obtained  to  within  14  percent  and 
C. . to  within  55  percent.  Finally,  the  moment  of  inertia 
coefficient  C.  is  obtained  to  within  0.02  percent. 


The  algorithm  is  sensitive  to  modeling  errors.  This 
is  reflected  in  the  measurement  residual  processes.  In  pro- 
cessing test  data,  the  residuals  can  indicate  that  the  correct 
projectile  model  is  not  being  used,  thus  helping  in  the  de- 
velopment of  a correct  model  for  a given  projectile.  The  EKF 
is  also  sensitive  to  the  initial  coefficient  errors  and  the 
initial  covariance  matrix  assumed  for  them.  If  these  errors 
and/or  their  covariance  are  too  large,  the  linearization 
assumed  by  the  filter  tray  be  inadequate  and  the  filter  may 
operate  poorly  or  diverge.  Fortunately,  engineering  judgment 
in  the  initialization  of  the  filter  has  been  sufficient  to 
insure  good  filter  performance  in  this  example. 

The  EKF  algorithm  requires  the  integration  of  the 
system  states  and  the  upper  triangle  of  a 29  > 29  covariance 
matrix.  This  is  done  by  trapezoidal  integration  at  a step 
size  of  0.05  millisecond.  Larger  step  sizes  lead  to  numerical 
instabilities  and  incorrect  covariance  propagation.  The  algo- 
rithm operates  in  approximately  300  thousand  bytes  of  memory 
on  an  IBM-360-1'15,  requiring  50  minutes  of  CPU  time  to  process 
the  50  sets  of  measurement  data.  An  operational  algorithm 
could  be  developed  which  would  improve  these  t imc  and  memory 
requirements  considerably  by  elimination  of  unnecessary  compu- 
tations and  storage  and  by  improvements  in  the  integration 
technique  and  program  organization. 


SECTION  V 


CONCLUSIONS 


The  principal  conclusions  of  this  report  are: 


The  feasibility  and  accuracy  of  the  extended  Kalman 
filter  for  ballistic  range  aerodynamic  coefficient 
extraction  has  been  demonstrated  for  a 30mm  projectile. 

For  the  example  used  in  this  work,  the  EKF  estimates  8 
of  the  projectile's  17  aerodynamic  coefficients  to 
within  10  percent  of  their  true  values.  The  remaining 
coefficients  have  such  a small  effect  on  the  trajectory 
that  they  cannot  be  estimated  this  accurately. 

Significant  coefficient  variations  with  mach  number 
should  be  modeled  and  can  be  estimated  by  the  EKF 
algorithm. 

The  EKF  accurately  estimates  the  rms  errors  associated 
with  each  parameter  estimate. 

Those  coefficients  which  cannot  be  estimated  from  the 
data  are  identified  by  the  EKF  algorithm,  and  the 
a priori  uncertainties  assumed  for  them  remain 
essentially  unchanged. 

The  EKF  technique  is  general  and  can  be  adapted  to 
various  situations  and  conditions.  It  incorporates  any 
a priori  knowledge  available  about  the  aerodynamic 
coefficients  which  might  be  available  from  design 
considerations  and  wind  tunnel  tests. 

The  absolute  parameter  estimation  accuracy  achievable 
is  sensitive  to  the  particular  trajectory  motion  that 
occurs,  the  measurement  noise  levels,  and  any  modeling 
errors  in  the  filter  design. 


SECTION  VI 


I * 


RECOMMENDATIONS 


r ■ 


The 

be  used  as  a 
evaluations , 
because 
d i t io  n s 


general  EKF  technique  demonstrated  herein  can 
tool  in  aerodynamic  model  development,  test 
and  ballistic  range  error-level  calibration 
of  its  sensitivity  to  modeling  errors  and  test  con- 
Other  possible  applications  include  the  reduction 


of  wind  tunnel  data  and  the  extraction  of  missile  and  air- 
craft aerodynamic  coefficients  from  onboard  and  external 
measurements . 


The  EKF  algorithm  should  be  compared  with  other 
current  data  reduction  techniques  such  as  the  Chapman-Kirk 
algorithm.  Development  of  an  operational  algorithm  would 
then  be  a natural  extension  of  the  current  work.  The  algo- 
rithm should  be  evaluated  for  projectiles  other  than  the 
30mm  round  used  in  these  investigations  to  assess  its  per- 
formance over  a range  of  aerodynamic  bodies.  Finally,  the 
expansion  of  the  current  model  to  the  asymmetric  airframe 
case  should  be  pursued. 
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APPENDIX  A 


COMPUTATION  OF  F(x(t),t) 


The  propagation  of  the  EKF  estimation  error  covariance 
matrix,  P(t),  between  measurements  requires  the  computation  of 
the  matrix  F(x(t),t)  as  indicated  in  Table  1 and  subsection  2.2 
This  matrix  contains  the  partial  derivatives  of  each  filter 
state  derivative  with  respect  to  each  filter  state,  evaluated 
at  the  current  state  estimate,  x ( t ) . This  appendix  gives  the 
equations  used  to  compute  this  matrix. 

PRELIMINARY  DEFINITIONS 

Table  A-l  summarizes  the  system  model  upon  which  the 
EKF  algorithm  design  is  based.  The  matrix  F(x(t),t)  is 
defined  with  respect  to  this  model  by 


F(x(t  ). t ) 


r af (x, t ) 

3x 


(A-l  ) 


Here,  the  element  in  the  ith  row  and  jth  column  of  F is 
given  by 


- 


3f i(x, t ) 


(A-2) 


where  each  function  f ^ is  defined  in  Tab  It  A-l,  with  its  cor- 
responding equation  given  in  Table  6.  The  expressions  giver, 
in  this  appendix  for  the  elements  of  F are  derived  by  direct 
differentiation  of  the  appropriate  functions.  The  only 
assumption  made  in  this  derivation  is  that  the  air  density, 

P,  is  independent  of  altitude,  z.  While  not  strictly  correct, 
the  error  induced  by  this  approximation  is  negligible. 

Intermediate  Constants 

The  following  variables  are  defined  for  convenience 
and  are  assumed  to  be  constant  for  the  purposes  of  computing 
th<  required  partial  derivatives. 

Tr  “ VTy 
r = d/2 


I am  = pA/2m 

ax  = pAd/2Ix 
ay  = pAd/2 I y 

Intermediate  Vari  abl  es  and  Their  Partial  Derivatives 

The  following  variables  (b^  through  and  through 

k^,, ) are  functions  of  the  system  states.  Their  definitions 

are  given,  followed  by  their  partial  derivatives  with  respect 
to  each  state  on  which  they  depend. 

1 J 

*■—  - TABLE  A-l  . EKF  STATE  VARIABLE  MODEL  SUMMARY 


Variable  bc 
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npu  npa  npu3 

( npa3b2u 

( np.x3b2v 

np  i3b2w 
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k2  = vVCNot 
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(v2+V2)CNa/V  + 

[V  ( 3v  + w )-v  (v  +w  )JCNa3/V 
vwC^/V  + vw(V2+u2)CNa3/V3 


’ariable  k. 


k3  = »’wCYpa 

k3u  = 3k3/3u  = pwb2uCYpa3 
k 3 v = 3k3/"V  = pwb2vrYp^3 


30 


= wVC 


Mu 


Mv 


4w 


Not 
3k  j/ 3u 

3kt/3v 

3 k j / 3 w 


UwCNa/V  " »w(v2^w2)CNfx3/V^ 
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k5 

= 

PV  CYpa 

k5u 

= 

3 k 5 / 3 u 

= pvh2uCYpu3 
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k12u 

3k12/Su  = 

PuClp/V 

k12v 

o k j £ / 9 v = 

PvClp/V 
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ELEMENTS  OF  F(x,t) 


This  subsection  gives  the  expressions  for  the  ele- 
ments o 1'  F(x,t)  in  terms  of  the  variables  defined  in  the  last 
section  and  in  subsection  2.2. 
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Row  1 of  F ( x , t ) 


Column , i 3f , / lx . 
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Pow  2 of  F ( x , t ) 
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Row  3 of  F ( x , t ) 
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Row  4 of  F(x , t ) 
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Rcnv  6 of  F ( x , t ) 


Column,  i 


3V3xi 


-am(k4u+rk5u)  + 0 

'am(k4v+rk5v')  + * sinf4 

-a  (k.  +rk,  ) 
m 4w  5w 


g sinfi  + v oos6 


v s i n 0 


-a  rkc 
ri  op 


-a  wV 
m 

-a  rpv 
m H 


-a  wVe 
m 

-a  rpve‘ 
m 


57 


* 

(• 


I 


Row  7 of  F ( \ . t ) 


• 

Co  1 urnn  , i 

of 7/?x. 

■ 

* \ 

1 

0 

r :■ 

9 

0 

r-  * 

1 h 

10 

1 

11 

0 

1 

29 

0 

« 
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Row  9 of  F(x , t ) 
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Row  12  of  F ( x , t ) 
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Rows  13  through  29  of  F(.\,t) 


Since  f.(A,t)  = 0 for  i 
in  rows  13  through  29  are  zero. 


.,29,  then  all  elements 


f 2 
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